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DISCONTINUOUS  DUAL-PRIMAL  MIXED  FINITE  ELEMENTS  FOR  ELLIPTIC 

PROBLEMS 

CARLO  L.  BOTTASSO*,  STEFANO  MICHELETTI*,  AND  RICCARDO  SACCO* 

Abstract.  We  propose  a  novel  discontinuous  mixed  finite  element  formulation  for  the  solution  of 
second-order  elliptic  problems.  Fully  discontinuous  piecewise  polynomial  finite  element  spaces  are  used 
for  the  trial  and  test  functions.  The  discontinuous  nature  of  the  test  functions  at  the  element  interfaces 
allows  to  introduce  new  boundary  unknowns  that,  on  the  one  hand  enforce  the  weak  continuity  of  the 
trial  functions,  and  on  the  other  avoid  the  need  to  define  a  priori  algorithmic  fluxes  as  in  standard 
discontinuous  Galerkin  methods.  Static  condensation  is  performed  at  the  element  level,  leading  to  a 
solution  procedure  based  on  the  sole  interface  unknowns. 

The  resulting  family  of  discontinuous  dual-primal  mixed  finite  element  methods  is  presented  in  the 
one  and  two-dimensional  cases.  In  the  one-dimensional  case,  we  show  the  equivalence  of  the  method 
with  implicit  Runge-Kutta  schemes  of  the  collocation  type  exhibiting  optimal  behavior.  Numerical 
experiments  in  one  and  two  dimensions  demonstrate  the  order  accuracy  of  the  new  method,  confirming 
the  results  of  the  analysis. 

Subject  classification.  Applied  and  Numerical  Mathematics 

Key  words,  finite  element  method,  mixed  methods,  discontinuous  Galerkin,  Petrov- Galerkin, 
elliptic  problem 

1.  Introduction  and  Motivation.  We  consider  the  following  classical  model  problem: 

f  -div(zA7u)  —  f  in  0  C  K2,  ^  ^ 

|  u  =  5  on 

where  H  is  a  bounded  polygon  with  Lipschitz  boundary  V  =  90, ,  and  v,  g  are  given  functions. 
Dirichlet  boundary  conditions  are  considered  only  for  the  sake  of  simplicity,  but  the  extension  to 
the  case  of  more  general  boundary  data  is  straightforward.  Definitions  of  functional  and  geometrical 
entities  are  given  in  section  2. 

In  view  of  the  approximation  of  (1.1),  we  let  {Th}h>o  be  a  family  of  triangulations  of  Q  and,  for 
any  7/i,  we  denote  by  T  any  element  thereof.  Introducing  the  auxiliary  unknown  a  =  problem 
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(1.1)  can  be  reformulated  in  weak  form  over  each  element  T  of  Th  as 


v  1a-rdx  +  /  udivrdx-  /  ur-nTds  =  0, 

h  Jot  (1.2) 

a  ■  V_v dx  —  /  va-nTds  —  /  fvdx, 

JdT  Jt 

where  r  and  v  are  smooth  vector  and  scalar  functions,  respectively,  on  T,  while  nT  is  the  unit  outward 
normal  vector  to  dT. 

A  unified  framework  for  the  analysis  of  the  Discontinuous  Galerkin  (DG)  method  applied  to  elliptic 
problems  in  the  form  (1.2)  was  provided  in  ref.  [1].  We  shall  follow  the  same  path  in  the  following, 
since  the  proposed  method  fits  well  in  that  framework.  The  DG  formulation  at  the  discrete  level  reads 
(see  also  [9]):  find  Uh  G  Uh  and  g_h  G  such  that  VT  G  % 


Uh  divr h  dx  —  y;  huTh  *  nT  ds  =  0  Vrft  G  S(T), 


eedT''?- 


/  a.h  *  Vvfc  dx-Y]  vhha  • nTds-  /  fvh 
Jt  -TZrJe  Jt 


dx 


Vvh  G  C/(T), 


e£dT  — 


(1.3) 


where  S(T)  and  U ( T )  are  two  sets  of  smooth  vector  and  scalar  functions,  typically  polynomials,  defined 
on  each  element  T,  while  and  Uh  are  the  corresponding  global  finite  element  spaces.  We  are  now 
dealing  with  a  Galerkin  method  where  the  values  of  the  unknown  fields  on  each  edge  e  of  <9T,  noted 
Uh\e_  and  vh\ei  have  been  replaced  by  the  numerical  fluxes  hu  and  ha.  In  practical  implementations, 
hu  and  ha  are  chosen  as  functions  of  the  internal  unknowns  Uh  and  ah  on  two  neighboring  triangles. 
In  particular,  it  can  be  shown  that  most  methods  reported  in  the  literature  can  be  classified  based  on 
the  particular  choice  of  the  numerical  fluxes  hu  and  ha ,  choice  often  inspired  by  ideas  borrowed  from 
finite  volume  methods.  The  interested  reader  can  consult  ref.  [1]  for  exact  details  on  the  expressions 
of  the  fluxes  in  the  different  cases. 

It  is  clear  that  any  particular  choice  on  the  nature  of  the  numerical  fluxes  will  affect  the  resulting 
scheme,  eventually  determining  its  numerical  characteristics,  stability,  accuracy  as  well  as  the  pattern 
of  the  stiffness  matrix.  This  approach  to  the  choice  of  the  numerical  fluxes  is  somewhat  unsatisfactory, 
in  the  sense  that  it  seems  to  be  a  rather  artificial  and  dogmatic  process,  in  some  sense  a  “violence” 
to  the  underlying  weak  form  (1.2).  In  this  work,  we  shall  therefore  depart  from  the  classical  approach 
of  defining  a  priori  the  expressions  of  the  numerical  fluxes.  Hence,  the  interface  fields  hu  and  h „  •  nT 
are  assumed  to  represent  new  additional  problem  unknowns ,  and  we  let  the  method  implicitly  define 
their  values.  In  doing  so,  we  obtain  a  genuine  Petrov-Galerkin  version  of  the  DG  method. 

This  formulation  generalizes  the  method  described  in  refs.  [2,  3,  4,  5]  to  multiple  spatial  dimen¬ 
sions.  In  the  one-dimensional  case,  the  method  can  be  shown  to  be  equivalent  to  collocation- type 
implicit  Runge-Kutta  (RK)  schemes.  In  particular,  when  used  in  conjunction  with  Gauss-Legendre 
quadrature,  the  finite  element  method  here  proposed  corresponds  to  the  Kuntzmann-Butcher  (Gauss) 
RK  scheme,  which  enjoys  well  known  properties  of  optimality .  This  provides  a  strong  motivation  for 
seeking  an  extension  of  those  ideas  to  the  multi-dimensional  case,  extension  that  is  attempted  in  this 
work  for  the  elliptic  boundary- value  problem. 

The  paper  is  laid  out  as  follows.  Section  2  presents  the  proposed  method  for  the  classical  problem 
of  the  Laplacian.  Section  3  deals  with  the  one- dimensional  formulation  of  the  present  method  and 
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with  its  connections  with  the  time  finite  element  and  RK  methods  of  refs.  [2,  3,  4,  5].  In  doing  so, 
we  shall  temporarily  abandon  the  simple  linear  model  problem,  and  present  the  method  for  a  general 
non-linear  ODE  problem.  Section  4  is  devoted  to  the  description  of  the  method  in  the  two-dimensional 
case.  Numerical  results  are  discussed  in  section  5,  where  the  proposed  method  is  validated  on  one  and 
two-dimensional  test  cases.  Finally,  the  conclusions  are  discussed  in  section  6. 

2.  The  Discontinuous  Dual-Primal  Mixed  Method  for  a  Model  Problem.  We  are  now 

ready  to  define  the  proposed  Discontinuous  Dual-Primal  Mixed  (DDPM)  discretization  of  the  one- 
element  weak  formulation  of  problem  (1.1):  find  Uh  G  g_h  G  £/*,  A^  G  A h,g  and  G  A h  such  that 
VT  e  Th  we  have 

/  v~l°h  'Thdx+  /  uh  di vrh  dx  -  Xhrh  •  nT  ds  =  0  Vr^  G  W(T), 

<  JT  Jt  jdT  (2.1) 

/  °h'Vvhdx-  /  SrVhfJ'h  ds  =  /  fvhdx  Vvh  €  V(T). 

JT  JdT  JT 

The  meaning  of  the  symbols  is  as  follows.  For  any  %  6  { Th}h>o ,  let  Eh  be  the  associated  set  of 
edges.  For  any  edge  e  of  Eh,  we  denote  by  ne  the  unit  normal  vector  to  e,  directed  according  to  a 
uniquely  fixed  arbitrary  convention.  Namely,  ne  coincides  with  the  unit  outward  normal  vector  on  the 
boundary  T  if  efl  T  ^  0,  otherwise  ne  coincides  with  one  of  the  two  outward  normals  to  the  triangles 
sharing  the  edge  e  if  eDT  =  0.  Then,  for  any  T  G  7*  we  denote  by  nT  the  unit  outward  normal  vector 
to  dT  and  define  the  sign  function  St  =  ne  •  nT  =  ±1.  Through  the  use  of  <St,  we  can  associate  the 
interface  degrees  of  freedom  with  mesh  edges,  so  that  the  sign  function  will  determine  whether  a  given 
flux  }ih  is  entering  of  exiting  from  a  given  triangle.  For  simplicity,  we  shall  omit  in  the  following  the 
dependence  of  n  on  the  edge  e  which  will  be  understood. 

Given  the  fact  that  the  edge  fluxes  are  now  unknowns  to  the  problem,  we  can  not  choose  test  and 
trial  functions  within  the  same  spaces,  as  for  the  DG  method  (1.3).  Therefore,  we  are  now  dealing 
with  a  Petrov- Galer kin  formulation.  The  global  finite  element  spaces  for  the  internal  variables  ah  and 
Uh  are  defined  as 


=  H  e  L2(Q)  I  vh\T  €  E(T)  VT  €  %} , 

Uh  =  {vh  €  L2(ft)  |  vh\T  €  U{T)VT  €  %} , 

where  L2(Q)  =  (L2(n))2.  Considering  now  the  interface  fields  A h  and  /x^,  we  note  that  they  represent 
the  traces  of  u  and  a  *  nT  on  dT ,  respectively.  Using  the  language  of  computational  mechanics,  the 
new  unknowns  act  as  “mixed  connectors”  between  neighboring  triangles  (see  [16]  and  [18],  Chapter 
13).  We  let  A(e)  indicate  a  set  of  scalar  smooth  functions  defined  on  each  edge  e.  Then,  we  define  the 
(global)  finite  dimensional  spaces 

A h  =  {vh  £  L2(Eh)  |  rjhle  £  A(e)  Ve  G  Eh}  , 

A  h,g  =  {r)h  ^Ah\r]h^  VgiienT^  0}, 

where  V  is  the  L2-projection  operator  on  A(e).  Finally,  we  introduce  the  local  and  global  finite  element 
spaces  for  the  test  functions  r  and  v.  Denoting  by  W_(T)  and  V  (T)  two  sets  of  smooth  vector  and 
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scalar  functions  defined  on  each  element  T,  we  have 

wh  =  {rh  €  I2(ft) \rh\T  €  W(T) VT  €  Th)  , 

Vi,  =  {w*  G  L2(Sl)  |  wA|T  e  V”(T)  VT  G  Tft}  . 

The  choice  of  the  functional  spaces  will  be  detailed  in  the  following  sections  for  the  one  and  the 
two-dimensional  case. 

We  note  that  the  method  is  completely  conservative  in  the  sense  of  ref.  [1],  since  the  interface 
unknowns  are  associated  with  edges  and  are  consequently  identical  for  two  triangles  sharing  the  edge. 
This  last  property  implies  the  following  conservation  statement 

V  Hhds- h  /  fdx  -  0, 

seouJs.  Jv 

where  U  is  the  union  of  some  collection  of  elements,  and  Vh  was  taken  to  be  identically  unity  in  the 
second  of  (2.1). 

The  DDPM  approach  can  be  viewed  as  a  consistent  hybridization  of  the  dual-primal  mixed  (DPM) 
method  introduced  in  ref.  [12].  Indeed,  the  spirit  of  the  standard  hybridization  procedure  (see  [7], 
Chapter  V)  is  to  relax  the  continuity  of  some  of  the  unknowns  (a  -n  in  the  dual  hybrid  philosophy  or 
u  in  the  primal  hybrid  philosophy)  and  then  to  enforce  it  back  through  the  introduction  of  suitable 
inter-element  Lagrange  multipliers.  Precisely,  u\e  represents  the  Lagrange  multiplier  in  the  case  of 
dual  hybrid  methods,  while  a  -  n \L  is  the  Lagrange  multiplier  in  the  case  of  primal  hybrid  methods. 
In  the  present  case,  we  relax  the  continuity  of  both  u  and  g_  •  n  in  the  interior  fields  and  let  the 
interface  variables  connect  neighboring  elements.  The  interface  variables  A  and  fi  do  not  play  the 
role  of  Lagrange  multipliers  in  our  case  but,  nevertheless,  enjoy  appealing  properties.  Physically,  they 
have  the  clear  meaning  of  u\e  and  cr  •  n|e,  respectively,  while  numerically  they  exhibit  a  higher  rate  of 
convergence  than  the  internal  variables,  as  it  will  be  shown  in  the  following.  This  is  also  a  feature  of 
the  Lagrange  multipliers  in  the  case  of  standard  hybrid  methods  ([7],  Chapter  V). 

Alternatively,  one  could  note  that,  by  gluing  the  test  functions  together  at  neighboring  triangles, 
each  internal  edge  receives  equal  and  opposite  contributions  from  the  two  elements  using  it,  and 
therefore  the  boundary  unknowns  are  eliminated.  This  way  one  recovers  from  (2.1)  the  DPM  method. 

3.  The  DDPM  Method  in  the  One-Dimensional  Case.  We  abandon  now  for  the  moment 
the  linear  model  problem  (1.1),  and  we  consider  a  generic  system  of  first  order  ODE’s 

£  =  lU,x),  (3.1) 

where  z,f_  G  R,J  and  D  is  the  number  of  state  variables.  Suitable  initial  or  two-point  boundary 
conditions  complement  problem  (3.1).  In  the  case  of  the  linear  model  problem,  we  clearly  have 
z  =  (ii,  a)T  and  /  =  (a,  — /)T.  However,  since  the  discussion  can  be  easily  completed  for  the  more 
general  case,  we  consider  (3.1)  instead  of  the  one-dimensional  version  of  (1.1)  in  the  remainder  of  this 
section. 

We  assume  0  =  (0, X)  and  let  0  =  <  X\  <  ...  <  xn~\  <  xn  =  X  be  a  given  partition  Th  of 

f }  into  n  >  1  intervals  T{  —  [xi,Xi+ 1]  of  size  hi,  i  —  0, . . .  ,n  —  1,  where  are  the  nodes.  Note 

that  in  the  one- dimensional  case  Sh  =  {xi}f-0  and  each  edge  e  degenerates  into  one  single  node. 
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The  statement  of  the  present  finite  element  method  is  as  follows:  find  (zh,  Xh)  £  (. Zh  x  Ah)  such 
that  for  each  T*  £  Tk->  i  =  0, . . .  ,  n  -  1,  the  following  holds 

/  zh‘w'hdx+  l(zh,x)-whdx-(Xw^wh(xi+1)~Xi-wh{xi))  =  0 

JTi  JTi 

VwheW(TiY  (3.2) 


For  any  T  G  7i,  we  define  the  local  finite  element  spaces  for  the  internal  fields  zh  and  for  the  test 
functions  as  follows: 


Z(T)  =  (¥k(T))D,  W(T)  =  (Pk+i(T))D, 
with  k  >  0.  The  local  finite  element  space  for  the  interface  unknowns  A^  is 

A (dTi)  —  {A^,  A^j},  i  =  0, . . .  ,  n  —  1, 

where  for  any  function  A  £  A(dT{)  we  let  A^  =  A  (a?*)  and  Ai+1  =  X(xi+i),  i  =  0, . . .  ,ra— 1.  Furthermore, 
we  set 


U{T)  =  U(T)  x  A  (AT),  V(T)  =  W(T). 

The  finite  element  space  is  then 

BPfc(T)  =  {(z,  Ajtfi)  |  U,  A;  w)  €  tf(T)  x  V(^)  VT  €  % }.  (3.3) 

The  global  finite  element  spaces  for  the  internal  fields  and  the  test  functions  are 

Zh  =  {zh  €  (L2(H))D  |  zh\T  e  Z(T)  vr  €  Tft}  , 

T£h  =  K  e  (i2(^))D  I  6  W.(T)YT  e  %}  , 

while  the  global  spaces  for  the  interface  unknowns  are 

Ah  =  {{Ai}?=o|Ai€RD,*  =  0, 

where  functions  are  defined  only  at  the  nodes  of  Sh •  Finally,  gathering  the  previous  definitions,  we 
obtain  the  global  finite  element  space  as 

mkh  =  (ZhxAh)xWh.  (3.4) 

It  is  clear  that  both  initial  and  boundary-value  problems  can  be  solved  within  this  framework. 
For  an  initial  value  problem,  one  can  solve  marching  sequentially  in  an  element-by-element  fashion,  as 
with  any  time  stepping  procedure.  We  have:  dim(Z(T))  =  D(h  + 1),  dim(W(T))  =  D(k  +  2)  VT  £  %- 
Therefore,  for  any  T  £  Th ,  (3.2)  leads  to  a  system  of  D(k  +  2)  equations  in  the  D(k  4-  1)  4-  2D 
local  unknowns  zh\T,  A/Jt?  leaving  space  for  the  D  initial  conditions  that  complement  problem  (3.1). 
A  boundary- value  problem  leads  to  a  global  discrete  problem  defined  over  0.  In  this  case  we  have 
nD(k+ 2)  equations  in  the  nD(k+l)~\~ 2nD  unknowns.  However,  considering  that  Ai+1  appears  in  both 
neighboring  elements  Ti  and  Tj+i ,  the  total  unknown  count  is  reduced  to  nD(k  + 1)  +  2 nD  —  (n  —  l)D, 
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which  leaves  once  again  space  for  D  two-point  boundary  values.  In  both  cases,  the  internal  field  zh 
can  be  eliminated  at  the  element  level  in  favor  of  the  interface  unknowns  \h ,  and  then  recovered  at 
convergence. 

We  note  that  the  scheme  presents  two  jump  discontinuities  for  each  finite  element,  namely  one  at 
X{  and  the  other  at  Xi+ 1,  since  A*  ^  zh(xi )  and  Ai+1  ^  zh(xi+ 1).  For  this  reason  it  was  termed  the 
Bi-Discontinuous  (BD)  scheme  in  ref.  [3].  It  is  clear  that  such  symmetric  treatment  of  the  element 
boundaries  implies  no  privilege  in  the  direction  of  the  flow  of  information,  and  it  is  therefore  best 
suited  for  two-point  boundary- value  problems. 

Quite  differently,  for  initial  value  problems  it  is  sometimes  desirable  to  have  stiffly  accurate 
schemes,  that  are  able  to  damp  the  higher  frequencies  components  from  the  computed  response  [11]. 
This  is  achieved  in  practice  by  forcing  a  lack  of  symmetry  in  the  scheme,  that  incorporates  the  knowl¬ 
edge  on  the  direction  of  flow  of  information  within  the  element.  It  is  possible  to  derive  such  schemes 
in  the  present  framework,  allowing  one  single  jump  discontinuity  at  Xi  while  enforcing  the  condition 
Ai+i  =  z.h(xi+i)-  One  obtains  in  this  case  the  so  called  Singly-Discontinuous  (SD)  schemes  [3]. 

3.1.  Equivalence  with  Runge-Kutta  Methods.  It  was  shown  in  ref.  [3]  and  then  again  with 
a  slightly  different  proof  in  ref.  [4],  that  the  finite  element  method  (3.2)  can  be  written  as  a  RK  process 
for  any  order  k.  RK  methods  are  probably  best  known  for  the  solution  of  initial  value  problems,  but 
clearly  can  also  be  used  for  the  solution  of  two-point  boundary-value  problems  by  assembling  the 
single  steps  to  yield  a  global  discrete  problem  defined  over  the  whole  computational  domain,  exactly 
as  in  the  finite  element  method.  Since  the  link  between  the  two  approaches  seems  to  be  useful  in  the 
characterization  of  the  proposed  method,  we  briefly  review  this  result  in  the  following. 

We  begin  by  selecting  a  quadrature  rule  for  the  evaluation  of  the  integrals  in  (3.2).  The  rule  is 
defined  by  s  abscissae  c*  and  weights  b{.  We  consider  the  case  s  =  k  - hi,  therefore  the  number  of 
quadrature  points  is  the  same  as  the  number  of  finite  element  nodes  of  the  trial  functions. 

Following  the  usual  FEM  practice,  one  expresses  the  discrete  equations  in  terms  of  the  nodal 
values.  Here  we  shall  take  a  different  approach,  and  assume  as  unknowns  the  values  of  the  finite 
element  trial  functions  at  the  quadrature  points .  This  is  the  key  idea  for  showing  the  link  existing 
between  finite  element  and  RK  methods.  The  change  of  unknowns  can  be  done  because  we  are  using 
a  quadrature  rule  that  uses  as  many  integration  points  as  the  number  of  finite  element  nodes,  as 
previously  said.  We  will  come  back  later  to  some  interesting  consequences  of  this  assumption. 

We  first  note  that  the  test  functions  can  be  expressed  as 

fc+2 

wh  =  ’YjTi~lai,  (3.5) 

i=l 


where  are  polynomials,  rGl,  and  ak  their  associated  amplitudes.  Next,  let  us  define  the  (s-bl)  x  s 
matrix  =  [t*-1  |r-^.] ,  with  i  =  1, . . .  ,  $  +  1,  j  =  1, . . .  ,  s,  the  £/s  being  s  real  values.  The  derivative 
of  with  respect  to  r  is  then  P £  =  [(i  —  l)r*"2|r-^] .  Note  that,  for  any  choice  of  the  £/s,  the  first 
row  of  P$  is  all  made  of  one’s,  and  the  first  row  of  P ^  is  all  made  of  zero’s.  For  =  Cj,  we  have  the 


,  with 


matrices  Pc  =  [cj  x] ,  P*c  =  [(i  -  1  )cj  2] .  We  also  define  the  s  x  s  diagonal  matrix  B  =  ^6*^ 
i  =  1, . . .  ,  5,  and  the  following  s-dimensional  vectors:  b  =  (6i , &2, . . .  >  bs)T,  c=  (c\ ,  c2, . . .  ,  cs)^,  and 


6 


0  =  (0,0,...  ,  Of,  1  =  (1, 1, ... ,  1)T,  li  =  (1, 0, . . . ,  of. 

We  note  that 


P'B  = 


0T 

Q' 


(3.6) 


Q  being  the  s  x  s  matrix  Q  =  [rl\ T=Cjbj]  =  \c)bj] .  The  derivative  of  Q  with  respect  to  r  is  given  by 

re  formula  is  of  order  s,  we  have  that 

Q'l  -  1.  (3.7) 


Q '  =  [i  clj  rbj] .  Furthermore,  if  the  quadrature  formula  is  of  order  s ,  we  have  that 


Equation  (3.7)  corresponds  to  the  so  called  “B  simplifying  assumption”  in  the  theory  of  RK  methods. 
The  discrete  weak  form  (3.2)  can  now  be  written  as 


((P'S)  ®  Id)  1 4-  hi((PcB)  ®  Id)  f_  =  (I  ®  Id)  Ai+1  -  (lx  ®  Id)  A*, 


(3.8) 


where  the  following  two  sD-dimensional  vectors  were  defined:  £  =  (i1,...  ,^5)r}  /  —  (/ (z1 , X{  + 
c\hi)i . .  .,[{zs,Xi  +  cshi))T.  Id  is  the  D  x  D  identity  matrix,  while  ®  denotes  the  tensor  Kronecker 
product  of  matrices,  so  that  (•)  ®  Id  ensures  that  all  degrees  of  freedom  of  the  vectorial  problem  are 
integrated  according  to  the  same  rule. 

Solving  for  z  and  Ai+1  and  using  (3.7),  we  have 


(i,Ai+if  = 


I  Mir  -  Q'~  Q) 

1  hibT 


(3.9) 


From  this  expression,  we  conclude  that  (3.2)  can  be  written  as  a  s-stage  RK  method  whose  tableau 
[10]  is  given  by 


1 bT  -Qf  Q 


(3.10) 


Finally,  we  test  the  Gauss-Legendre,  Lobatto  and  Radau  quadrature  rules  and  compute  the  cor¬ 
responding  tableaux  (3.10).  The  results  are  summarized  in  Table  (3.1). 

Table  3.1 

Correspondence  between  quadrature  rules  for  BD  finite  elements  and  RK  methods. 


Quadrature  Rule 

RK  Method 

Gauss 

Kuntzmann-Butcher 

Lobatto 

Lobatto  IIIB 

Radau-Left 

Radau  IA 

Table  (3.1)  shows  that  the  family  of  schemes  deriving  from  the  one-dimensional  finite  element 
method  considered  in  this  work  corresponds  to  some  well  known  RK  algorithms.  The  two  approaches 
differ  on  the  choice  of  the  unknowns:  the  finite  element  method  solves  for  the  nodal  values ,  while 
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the  RK  method  solves  for  the  quadrature  point  values.  However,  these  two  sets  of  values  are  simply 
related  by  a  linear  transformation.  Furthermore,  note  that  all  RK  schemes  appearing  in  the  table 
are  of  the  collocation  type,  where  polynomials  interpolate  the  two  interface  values  Ai5  Ai+i  and  the 
internal  stages  z  [10].  This  link  between  finite  elements  and  RK  processes  allows  a  somewhat  unique 
way  of  looking  at  the  discretization  procedure:  from  the  finite  element  point  of  view  we  see  jumps 
in  the  field  variables,  since  we  have  discontinuous  interpolations  of  the  internal  fields,  glued  together 
from  one  element  to  its  neighbor  through  the  presence  of  the  interface  unknowns.  Quite  differently, 
from  the  RK  point  of  view  we  have  continuous  interpolations  of  the  interface  unknowns  with  internal 
variables  at  the  collocation  points,  with  apparently  no  jump  discontinuities  in  the  solution.  It  should 
also  be  remarked  that  for  these  RK  methods,  maximal  order  is  in  general  attained  only  at  the  interface 
values,  and  not  at  the  internal  stages.  We  then  expect  to  observe  higher  rates  of  convergence  of  the 
interface  fields  also  in  the  multi-dimensional  generalization  of  the  method. 

The  same  table  shows  that  the  present  finite  element  formulation,  when  used  in  conjunction 
with  Gauss-Legendre  quadrature,  yields  the  Kuntzmann-Butcher  (or  Gauss)  RK  methods.  These  RK 
schemes  are  probably  the  ones  that  enjoy  the  greatest  number  of  numerical  properties,  and  are  in 
this  sense  optimal.  In  fact  they  are  of  maximal  order  (2s),  they  are  algebraically  stable,  and  also 
symplectic  when  applied  to  Hamiltonian  systems  [17].  Simplecticity  is  the  distinguishing  property  of 
Hamiltonian  systems,  and  its  numerical  preservation  has  a  strong  influence  on  the  behavior  of  the 
integration  scheme. 

As  a  final  point,  we  note  that  these  results  are  solely  based  on  having  assumed  a  quadrature  rule 
that  uses  as  many  points  as  the  number  of  finite  element  nodes.  In  reality,  working  from  the  point 
of  view  of  the  finite  element  method,  one  could  choose  to  use  a  greater  number  of  quadrature  points, 
for  example  in  the  hope  of  improving  the  accuracy  in  the  evaluation  of  strongly  non-linear  functions 
in  the  term  J  f(zh,x)  *  whdx .  It  can  be  shown  however,  that  increasing  the  number  of  quadrature 

points  can  in  reality  degrade  the  performance  of  the  method.  For  example,  using  a  larger  number  of 
points  destroys  the  symplectic  nature  of  the  scheme  [3]. 

4.  The  Two-Dimensional  Case.  For  two  spatial  dimensions,  the  choice  of  the  finite  element 
spaces  for  the  internal  and  interface  unknowns  is  a  straightforward  extension  of  those  introduced  in 
the  one-dimensional  case.  Precisely,  for  k  >  0  we  set 

W)  =  (F*(T))2 ,  U(T)  =  Pfc(T)  VT  €  %, 

A(e)=P*(e)  Ve€fh. 

The  choice  of  the  finite  element  test  spaces  is  less  obvious  and  requires  some  care.  Here  we 
focus  on  the  description  of  the  method  in  the  cases  k  —  0  and  h  —  1.  We  start  by  introducing  the 
Raviart-Thomas  finite  element  space  of  degree  k  [14] 

KT *(T)  =  (P k(T))2  ©  aP*(T)  VT  6  %. 

In  the  case  of  the  BP°  finite  elements,  we  set  W(T)  =  KT0(T)  and  V(T)  =  Pj(T).  The  global  finite 
element  spaces  W_h  and  Vh  are  then  simply  defined  as  the  products  of  the  corresponding  local  spaces. 
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Notice  that  W_h  and  Vh  are  the  discontinuous  counterparts  of  the  corresponding  test  spaces  that  have 
been  used  in  the  dual-primal  formulation  of  ref.  [12]. 

In  the  case  of  the  BP^  finite  elements,  we  set  W(T)  =  MTi(T)  and  V(T)  =  P2 (T)  0  B$(T)  for 
each  T  G  B^(T)  is  the  following  cubic  bubble  [15] 

£3  =  (<Po  ~  <Pi)(Pi  ”  <P2)((p2  -  <po), 

where  tpi  =  ipi(x ),  i  =  0, 1,2,  are  the  centroidal  coordinates  of  a  point  x  G  R2  with  respect  to  the 

vertices  of  T.  Enrichment  of  the  scalar  test  space  is  a  standard  procedure  in  hybrid  methods  when 

the  degree  of  the  polynomials  in  V(T)  is  even  (see  refs.  [15,  8]).  Actually,  one  can  check  that  without 

adding  the  cubic  bubble  the  rectangular  matrix  arising  from  the  integral  /  SrVh^h  ds  becomes  rank 

JdT 

deficient. 

As  far  as  the  implementation  is  concerned,  we  point  out  that  for  both  the  k  ~  0  and  k  =  1  cases 
it  is  possible  to  perform  a  static  condensation  of  the  internal  variables  in  favor  of  the  edge  unknowns, 
obtaining  a  linear  system  in  the  sole  interface  variables  A h  and  fi 

The  generalization  of  the  DDPM  method  to  higher  orders  will  be  more  extensively  addressed  in  a 
forthcoming  paper.  Let  us  just  mention  that  a  possible  general  recipe  for  constructing  the  PP£  finite 
element  spaces  for  even  k  can  be  given  as 

W(T)  =  BDFMk+\(T),  V(T)  =  P*+i(T)  VT  G  Th) 

where  BDFM^i  is  the  triangular  analogue  of  the  reduced  element  introduced  in  ref.  [6]. 

5.  Numerical  Results.  In  this  section  we  demonstrate  the  BP^  finite  elements  for  k  —  0  and 
k  =  1  through  the  solution  of  test  cases  in  both  one  and  two  spatial  dimensions. 

5.1.  One-Dimensional  Examples.  We  consider  the  following  two-point  boundary- value  model 
problem 


f  -(v(x)v!(x)y  =/(x),  0  <  £  <  1,  /51x 

\  u( 0)  -  u(  1)  =  0. 

In  all  numerical  experiments,  successively  finer  uniform  grids  of  size  h  =  2-J  ,  j  G  {1, . . .  ,9}  are  used. 
We  set  eu  =  u  —  Uh  and  ea  —  a  -  07*;  moreover  and  e ^  denote  the  functions  defined  only  at  the 
nodes  of  £h  such  that  e\ (x{)  =  u(xi)  -  A*  and  e^(Xi)  =  a(xi)  -  fiiy  i  —  0, . . .  ,  n. 

For  any  function  v  G  L2(H)  and  sufficiently  smooth  on  each  T  G  7*,  we  denote  by  the 

approximation  of  fT  v  dx  using  the  trapezoidal  rule.  Moreover,  for  each  T  G  Th  we  denote  by  {x* iT}*_0 
the  k  +  1  Gauss-Legendre  points  on  T.  In  order  to  measure  the  convergence  rate  of  the  BP  h  finite 
elements,  we  use  the  following  three  norms: 


if  k  =  0, 
if  k  =  1, 
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where  w  is  any  bounded  function  defined  at  the  Gauss  points  of  Th  and  rj  is  any  bounded  function 
defined  at  the  nodes  of  Th-  Clearly,  ||  *  11^0,0  is  a  discrete  L'2(Q.)  norm,  while  the  other  two  norms  are 
discrete  maximum  norms  over  the  set  of  Gauss  points  and  £&,  respectively.  The  quadrature  formula 
used  in  the  discrete  norm  is  accurate  enough  not  to  pollute  the  computed  accuracy,  and  at  the  same 
time  avoids  to  sample  the  solution  at  superconvergent  points.  Moreover,  observe  that  any  integral 
appearing  in  the  DDPM  discretization  of  (5.1)  has  been  computed  using  a  one-point  Gauss-Legendre 
quadrature  rule  (with  node  x [}  on  each  T  E  Th)  for  k  =  0  and  a  two-point  Gauss-Legendre  quadrature 
rule  (with  nodes  {xj>T}j=0  on  each  T  E  %)  for  k  =  1. 

The  following  symbols  are  used  in  the  graphs:  for  any  h1  symbol  refers  to  He^U^o^  (or 
||e<r|ko,n)5  while  symbols  and  ‘o’  denote  \\eu\\h^p,oo  (or  \\eff\\h,Gp, oo)  and  ||eA|U,oo  (or  ||e„||fc>0o), 
respectively. 

Test  Case  1  We  set  v(x)  =  1  and  f(x)  =  sin(7nc),  x  E  [0,1].  The  exact  solution  is  the  smooth 
function  u(x)  =  (1  /tt2)  sin  (to),  which  implies  cr(x)  =  u((x)  =  (l/7r)  cos(7rx).  We  show  in  figure  5.1 
the  error  curves  for  eu  and  (left)  and  ea  and  (right)  in  the  case  k  =  0.  The  corresponding  error 
curves  in  the  case  k  —  1  are  shown  in  figure  5.2. 


Fig.  5.1.  One  spatial  dimensions,  test  case  1.  Approximation  errors  eu  and  e\  (left)  and  ea  and  e ^  (right)  for 
DSP  k  finite  elements. 

The  plots  show  second-order  convergence  for  the  interface  unknowns  in  the  case  k  =  0,  while  the 
internal  unknowns  are  only  first-order  accurate.  Superconvergence  of  the  internal  variables  can  be 
observed  at  the  Gauss  point.  For  the  DSP^  finite  elements,  the  error  curves  show  that  the  interface 
unknowns  are  fourth-order  convergent,  the  internal  fields  are  second-order  accurate  while  a  third-order 
convergence  rate  is  exhibited  by  the  internal  unknowns  at  the  two  Gauss  points  of  each  element.  All 
these  results  are  in  agreement  with  the  theory  of  RK  methods. 


10 


Fig.  5.2.  One  spatial  dimensions,  test  case  1.  Approximation  errors  eu  and  e\  (left)  and  ea  and  e ^  (right)  for 
OP h  finite  elements. 

Test  Case  2  This  test  problem  demonstrates  the  performances  of  the  DDPM  method  in  the  presence 
of  strong  solution  gradients.  The  coefficient  v  and  the  source  term  /  are  chosen  as  in  ref.  [13]  (see 
also  ref.  [8],  section  4.4.1,  for  further  comments),  and  read 

u(x)  =  —  4-  a(x  —  x )2, 
a 

/(x)  —  2  (l  4-  a(x  -  x)  {tan_1(a(x  —  of))  +  tan-1  (ax)})  . 

The  exact  solution  is 

u(x)  =  (1  -  x)  (tan“1(a(x  -  x))  +  tan-1  (ax))  , 

°(x)  =  (^  +a:(a:-x)2)  [-  tan_1(a(x  -  x))  -  tan-^aar)  +  (1  ~  x)  ’ 

where  a(x)  =  i/(x)u/(x).  Furthermore,  we  set  a  =  100  and  x  =  0.36388. 

We  show  in  figure  5.3  the  error  curves  for  eu  and  e\  (left)  and  ea  and  (right)  in  the  case  k  =  0. 
The  corresponding  error  curves  in  the  case  k  =  1  are  shown  in  figure  5.4. 

The  numerical  results  clearly  show  the  increased  “roughness”  of  the  problem  when  compared 
with  the  corresponding  error  curves  for  test  case  1.  Nevertheless,  all  variables  attain  the  expected 
convergence  rate  as  h  — y  0.  Figure  5.5  shows  the  exact  solution  u  (solid  line)  plotted  against  the  finite 
element  solution  A h  computed  with  BP^  elements  and  h  =  1/20.  Notice  the  sharp  internal  layer, 
resolved  within  one  element. 

5.2,  Two-Dimensional  Examples.  For  the  two-dimensional  case,  we  consider  (1.1)  in  the 
square  domain  =  (— 1,1)2  with  g  =  0.  The  exact  solution  is  in  this  cas eu(x,y)  =  (x2-l)(?y2-l)  and 
t z(x,y )  =  2  ( x(y 2  —  l),y(x2  —  1))T;  the  right-hand  side  /  is  computed  accordingly.  In  all  numerical 
experiments,  we  use  a  uniform  grid  of  isosceles  right  triangles  of  size  hx  =  hy  =  h ,  where  h  takes  on 
the  values  21-^/5,  j  e  {0,1, 2, 3}  for  the  BP^  case  and  {2/5,2/10, 2/20, 2/30}  for  the  BP^  case. 
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Fig.  5.3.  One  spatial  dimensions ,  test  case  2.  Approximation  errors  eu  and  e\  (left)  and  e a  and  e ^  (right)  for 
DP  h  finite  elements. 


Fig.  5.4.  One  spatial  dimensions ,  test  case  2.  Approximation  errors  eu  and  e\  (left)  and  ea  and  e ^  (right)  for 
DP  h  finite  elements. 


As  in  the  one- dimensional  examples,  we  denote  by  eu  —  u  ~uu  and  —  Q_  —  g_h  the  discretization 
errors  associated  with  the  internal  fields,  while  e\  —  V\  —  A h  and  =  Vfi  —  Hh  are  the  discretization 
errors  associated  with  the  interface  fields.  V  is  the  L2 -projection  operator  on  A(e),  for  each  e  €  £h- 
Furthermore,  for  any  T  G  Th  we  denote  by  |T|  the  measure  of  T,  by  ar,  r  —  0, 1,2  the  vertices  of  T, 
by  xc  T  the  centroid  of  T  and  by  —  0, 1, 2  the  three  Gauss  points  of  T  obtained  from  the  area 

coordinates  {2/3, 1/6, 1/6}  by  permutation.  In  order  to  measure  the  convergence  rate  of  the  JW h 
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Fig.  5.5.  One  spatial  dimensions,  test  case  2.  Exact  solution  (solid  line)  and  numerical  solution  X h  (‘o’)  computed 
using  ©IP h  finite  elements  (h  =  1/20/. 


if  k  =  0, 
if  k  =  1, 


where  v  is  any  sufficiently  smooth  function  on  T  belonging  to  L2(Q ),  r)  is  any  function  belonging  to 
L2(e)  for  each  ee  Eh  and  ||7/||o,e  is  the  L2-norm  on  edge  e.  Notice  that  ||  *  \\h,o,n  is  a  discrete  L2(Q ) 
norm,  ||  ■  ||/i,gp,oo  is  a  discrete  L°°(tt)  norm,  while  ||  *  \\h,-i/2,eh  is  a  discrete  H~1^2(Eh)  norm  (see  also 
ref.  [7],  Sect.  V.3).  For  vector  functions  v  —  {vi,V2)t •>  the  above  norms  are  defined 

/  \  1/2 
IHk,GP,oo  -  maxie{1>2}  ||vi|U, gp,oo. 

In  the  following  figures,  symbol  is  used  to  denote  ||eu||k)0>n  (or  UeJI/^o.n),  while  symbols  ‘ET 
and  ‘o’  denote  \\eu\\h)GP)00  (or  ||e^ |U,gp,oo)  and  |teA|U ,~i/2,£h  (or  ||eJA  -i/2,£h),  respectively. 

We  show  in  figure  5.6  the  error  curves  for  eu  and  e\  (left)  and  and  e ^  (right)  in  their  respective 
norms  for  the  case  k  =  0.  The  corresponding  error  curves  in  the  case  k  —  1  are  shown  in  figure  5.7. 

For  the  u  field,  the  plots  show  second-order  convergence  for  the  interface  unknowns  in  the  case 
k  =  0,  while  the  internal  unknowns  are  only  first-order  accurate.  Superconvergence  of  the  internal 
variables  can  be  observed  at  the  Gauss  point.  This  behavior  is  similar  to  the  one  observed  in  the 
one-dimensional  case.  For  the  a  field,  the  plots  show  once  again  second-order  convergence  for  the 


finite  elements,  we  introduce  the  following  discrete  norms: 


ITI  2  '  1/2 

=  (  Z  VZ/^r))2 


IM[/i,GP,oo  = 


T€Th  r= 0 

J  max|«(xCiT)|, 


max  max 
TeThje{ 0,1,2}' 

/  \  W 


IMk-i/2A=  ElsIlMlL 

\e££h  j 
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Fig.  5.6.  Two  spatial  dimensions.  Approximation  errors  eu  and  e\  (left)  and  e^  and  e ^  (right)  for  DP h  finite 
elements. 


Fig.  5.7.  Two  spatial  dimensions.  Approximation  errors  eu  and  e\  (left)  and  e^  and  (right)  for  DP h  finite 

elements. 


interface  unknowns,  while  the  internal  unknowns  are  only  first-order  accurate.  However,  for  this  field 
superconvergence  at  the  Gauss  point  is  not  observed. 

In  the  case  of  DP^  finite  elements,  the  error  curves  for  the  u  field  show  third  order  accuracy 
for  the  interface  unknowns,  and  second  order  accuracy  for  the  internal  unknowns.  For  the  g_  field, 
the  interface  unknowns  achieve  an  order  equal  to  approximately  5/2,  while  the  internal  variables  are 
second-order  accurate.  Once  again,  superconvergence  is  not  observed  at  the  Gauss  points. 
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In  all  cases,  the  interface  unknowns  are  consistently  more  accurate  than  the  internal  fields,  as 
expected  from  the  one- dimensional  analysis,  although  the  exact  order  achieved  by  the  various  fields 
does  not  seem  to  obey  an  obvious  rule.  As  far  as  the  order  of  convergence  of  A h  is  concerned,  the 
results  agree  with  the  classical  convergence  analysis  of  dual-hybrid  methods.  Actually  in  this  case 
one  obtains  0(hk+2)  when  the  Raviart-Thomas  finite  elements  of  degree  k  are  employed  (see  [7],  sect. 
V.3,  formula  (3.19)). 

Let  us  conclude  by  comparing  the  performance  of  the  DDPM  approach  to  the  performance  of  the 
standard  dual  method.  In  particular,  we  consider  the  B ¥°h  method  with  static  condensation  versus  the 
dual  method  using  lowest  order  Raviart-Thomas  finite  elements  (B>J  with  no  hybridization.  Table  5.1 
shows  the  results  obtained  for  the  example  here  considered.  The  meaning  of  the  various  entries  is  the 
following:  h  is  the  value  of  the  mesh  size  on  the  boundary  of  the  domain,  dim  is  the  dimension  of  the 
matrix  associated  with  the  linear  system,  nnz  is  the  number  of  non-zero  elements  of  the  same  matrix 
and  flops  is  the  number  of  floating-point  operations  as  computed  by  Matlab  using  sparse  operations. 
The  table  on  the  left  refers  to  the  BP°  case,  while  the  table  on  the  right  collects  the  results  for  the 
Bh  case.  In  the  first  case,  flops  includes  assembly,  solution  of  the  linear  system  and  recovery  of  the 
internal  fields,  while  in  the  second  case  it  considers  the  sole  assembly  and  solution  phases,  since  there 
are  no  recoveries  to  perform. 


Table  5.1 

Flop  counts  for  the  BP  h  and  B^  methods. 


h 

dim 

nnz 

flops 

h 

dim 

nnz 

flops 

2/5 

150 

662 

41424 

2/5 

135 

617 

26888 

2/10 

600 

2713 

224775 

2/10 

520 

2272 

209964 

2/20 

2400 

11060 

1835876 

2/20 

2040 

9752 

2993179 

2/40 

9600 

44267 

16867719 

2/40 

8080 

36112 

22407438 

The  tables  show  that  the  BP°  method  performs  better  than  the  B^  method  in  terms  of  flop-count 
as  the  dimension  of  the  linear  system  grows,  despite  the  fact  that  the  dimension  of  the  matrix  and 
the  number  of  its  non-zero  entries  are  greater  than  in  the  B^  case. 

6.  Conclusions.  We  have  presented  a  novel  mixed  dual-primal  finite  element  formulation  for 
elliptic  boundary- value  problems.  This  work  is  a  first  attempt  at  generalizing  discontinuous  finite 
element  formulations  for  generic  ODE  problems  that  were  developed  in  previous  papers.  The  one¬ 
dimensional  formulation  can  be  shown  to  lead  to  optimal  RK  methods,  and  therefore  provides  a  strong 
incentive  towards  the  generalization  to  multiple  space  dimensions. 

The  DDPM  method  uses  interpolations  of  the  unknown  fields  which  are  discontinuous  between 
neighboring  triangles,  in  the  same  spirit  as  the  DG  method.  However,  no  particular  expression  of  the 
numerical  fluxes  is  selected  a  priori  in  this  case.  Additional  unknowns  are  introduced  at  the  element 
interfaces,  that  act  as  connectors  gluing  neighboring  triangles  together.  The  discontinuous  internal 
unknowns  can  be  eliminated  at  the  element  level,  leaving  a  solution  scheme  in  the  sole  interface 
variables.  A  higher  rate  of  convergence  is  observed  for  the  interface  unknowns  with  respect  to  the 
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internal  variables,  as  expected  from  the  analysis  in  the  one-dimensional  case.  Numerical  examples 
have  been  presented  to  support  the  analysis  and  provide  numerical  evidence  of  the  characteristics  of 
the  method  here  investigated. 

Future  efforts  will  concentrate  on  the  extension  of  these  ideas  to  other  model  problems,  and  on 
comparisons  of  this  class  of  methods  with  alternative,  well  established  solution  procedures. 
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